home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / language / parallax / more_exa.tar / more / fft.p next >
Text File  |  1991-11-26  |  7KB  |  264 lines
  1. SYSTEM  fastfouriertrans;
  2.  
  3. (* ======================================================= *)
  4. (* Program for parallel Fast-Fourier transformation        *)
  5. (* ======================================================= *)
  6.  
  7.  
  8. CONST n  = 4;   (* number of inputvalues (power of 2)      *)
  9.  
  10. TYPE  komplex = RECORD
  11.                   r : REAL; (* real part *)
  12.                   i : REAL; (* imaginary part *)
  13.                 END;
  14.       vek1    = ARRAY [0..n-1] OF komplex;
  15.  
  16.  
  17. (* ------------------------------------------------------- *)
  18. (* connection of single processing elements                *)
  19. (*   the connection structure is a two-dimensioned array   *)
  20. (*   (like a grid with connections in all four directions) *)
  21. (* ------------------------------------------------------- *)
  22.  
  23.  
  24. CONFIGURATION  liste [n];
  25. CONNECTION     left:  liste[i]  ->  liste[i-1].right;
  26.                right: liste[i]  ->  liste[i+1].left;
  27.  
  28.  
  29.  
  30. (* ------------------------------------------------------- *)
  31. (* variables :   (for main-program)                        *)
  32. (* ------------------------------------------------------- *)
  33.  
  34. SCALAR j, k, logn, h, p, q            : INTEGER;
  35.        a, b                           : vek1;
  36.        tmp2                           : komplex;
  37. VECTOR c, c_plus_p, z, 
  38.        tmp1, tmp5, 
  39.        tmp6, tmp7                     : komplex;
  40.        tmp3                           : INTEGER;
  41.  
  42.  
  43. (* ------------------------------------------------------- *)
  44. (* The procedure multplex computes the multiplication of   *)
  45. (* two complex numbers                                     *)
  46. (* ------------------------------------------------------- *)
  47.  
  48. PROCEDURE multkomplex(VECTOR a,b : komplex;VECTOR VAR c : komplex);
  49.  
  50. BEGIN
  51.   c.r := (a.r * b.r) - (a.i * b.i);
  52.   c.i := (a.r * b.i) - (a.i * b.r);
  53. END multkomplex;
  54.  
  55.  
  56.  
  57. (* ------------------------------------------------------- *)
  58. (* The procedure addkomplex computes the sumary of two     *)
  59. (* complex numbers                                         *)
  60. (* ------------------------------------------------------- *)
  61.  
  62. PROCEDURE addkomplex(VECTOR a,b : komplex;VECTOR VAR c : komplex);
  63.  
  64. BEGIN
  65.   c.r := (a.r + b.r);
  66.   c.i := (a.i + b.i);
  67. END addkomplex;
  68.  
  69.  
  70.  
  71. (* ------------------------------------------------------- *)
  72. (* The procedure subkomplex computes the subtraction of two*)
  73. (* complex numbers                                         *)
  74. (* ------------------------------------------------------- *)
  75.  
  76. PROCEDURE subkomplex(VECTOR a,b : komplex;VECTOR VAR c : komplex);
  77.  
  78. BEGIN
  79.   c.r := (a.r - b.r);
  80.   c.i := (a.i - b.i);
  81. END subkomplex;
  82.  
  83.  
  84.  
  85.  
  86. (* ------------------------------------------------------- *)
  87. (* The procedure multkomplexspez computes the power of a   *)
  88. (* whole number with a complex number                      *)
  89. (* ------------------------------------------------------- *)
  90.  
  91. PROCEDURE multkomplexspez(VECTOR a: komplex; VECTOR b: INTEGER;
  92.                           VECTOR VAR c : komplex);
  93.  
  94. VECTOR i    : INTEGER;
  95.        tmp : komplex;
  96.  
  97. BEGIN
  98.   IF b = 0 THEN
  99.     c.r := 1.0;
  100.     c.i := 0.0;
  101.   ELSE
  102.     tmp := a;
  103.     FOR i := 2 TO b DO
  104.       multkomplex(a,tmp,a);
  105.     END;
  106.     c := a;
  107.   END;
  108. END multkomplexspez;
  109.  
  110.  
  111.  
  112.  
  113.  
  114. (* ------------------------------------------------------- *)
  115. (* The procedure w_power computes the n-th unionradix w    *)
  116. (* with w = e^((2*PI)/n)                                   *)
  117. (* ------------------------------------------------------- *)
  118.  
  119. PROCEDURE w_hoch (SCALAR p : INTEGER; SCALAR VAR x : komplex);
  120.  
  121. SCALAR   tmp       : REAL;
  122.  
  123. BEGIN
  124.   tmp    := (2.0 * PI * FLOAT(p))/FLOAT(n);
  125.   x.r := Cos(tmp);
  126.   x.i := Sin(tmp);
  127. END w_hoch;
  128.  
  129.  
  130. (* ------------------------------------------------------- *)
  131. (* The procedure Reverse transforms a number :             *)
  132. (*   the transfering number would be decomposed in its     *)
  133. (*   binary representation, turned back and transfering    *)
  134. (*   back as decimal number.                               *)
  135. (* Important: representation in log(n)-Bit                 *)
  136. (* ------------------------------------------------------- *)
  137.  
  138. PROCEDURE reverse (SCALAR k : INTEGER; SCALAR VAR erg : INTEGER);
  139.  
  140. SCALAR tmp3             : REAL;
  141.        tmp1,tmp2        : INTEGER;
  142.  
  143. BEGIN
  144.   tmp1 := 0;
  145.   tmp3 := 0.0;
  146.   WHILE tmp1 <> (logn) DO
  147.     tmp2 := k DIV 2**(logn-1-tmp1);
  148.     IF tmp2 = 1 THEN
  149.       tmp3 := tmp3 + FLOAT(2**(tmp1));
  150.       k := k - 2**(logn-1-tmp1);
  151.     END;
  152.     tmp1 := tmp1 + 1;
  153.   END;
  154.   erg := TRUNC(tmp3);
  155. END reverse;
  156.  
  157.  
  158. PROCEDURE reverse2 (VECTOR k : INTEGER; VECTOR VAR erg : INTEGER);
  159.  
  160. VECTOR tmp3                  : REAL;
  161.        tmp1,tmp2             : INTEGER;
  162.  
  163. BEGIN
  164.   tmp1 := 0;
  165.   tmp3 := 0.0;
  166.   WHILE tmp1 <> (logn) DO
  167.     tmp2 := k DIV 2**(logn-1-tmp1);
  168.     IF tmp2 = 1 THEN
  169.       tmp3 := tmp3 + FLOAT(2**(tmp1));
  170.       k := k - 2**(logn-1-tmp1);
  171.     END;
  172.     tmp1 := tmp1 + 1;
  173.   END;
  174.   erg := TRUNC(tmp3);
  175. END reverse2;
  176.  
  177.  
  178.  
  179. (* ------------------------------------------------------- *)
  180. (* The procedure Input read input vector a                 *)
  181. (* ------------------------------------------------------- *)
  182.  
  183. PROCEDURE Input (SCALAR VAR a : vek1);
  184.  
  185. SCALAR i  : INTEGER;
  186.  
  187. BEGIN
  188.     WriteString(' Enter Input-Vector a :');
  189.     WriteLn;
  190.     FOR i := 0 TO n-1 DO
  191.       WriteString(' a['); WriteInt(i,1);
  192.       WriteString('] : ');
  193.       ReadReal(a[i].r);
  194.       a[i].i := 0.0;
  195.     END;
  196.     WriteLn;
  197. END Input;
  198.  
  199.  
  200.  
  201.  
  202.  
  203. (* ------------------------------------------------------- *)
  204. (* The procedure output displays solution-vector b         *)
  205. (* ------------------------------------------------------- *)
  206.  
  207. PROCEDURE output (SCALAR b : vek1);
  208.  
  209. SCALAR i,i2  : INTEGER;
  210.  
  211. BEGIN
  212.     WriteString(' Solution : ');
  213.     FOR i := 0 TO n-1 DO
  214.       reverse(i,i2);
  215.       WriteLn;
  216.       WriteFixPt(b[i2].r,7,2);
  217.       WriteFixPt(b[i2].i,7,2);
  218.     END;
  219.     WriteLn;
  220. END output;
  221.  
  222.  
  223.  
  224.  
  225. (* ------------------------------------------------------- *)
  226. (* start of main-program                                   *)
  227. (* ------------------------------------------------------- *)
  228.  
  229. BEGIN
  230.   Input(a);
  231.   LOAD (c,a);
  232.   logn := TRUNC (Ln(FLOAT(n))/Ln(2.0));
  233.   FOR h := logn-1 TO 0 BY -1 DO
  234.     p := 2**h;
  235.     q := TRUNC(FLOAT(n)/FLOAT(p));
  236.     PARALLEL 
  237.       w_hoch(p,tmp2);
  238.       z := tmp2;
  239.       tmp1 := c;
  240.       reverse2 (DIM1,tmp3);
  241.       multkomplexspez (z,(tmp3 MOD q),tmp5);
  242.       c_plus_p.r := c.r;
  243.       c_plus_p.i := c.i;  
  244.       tmp7.r := 0.0;
  245.       tmp7.i := 0.0; 
  246.       PROPAGATE.left^p(c_plus_p.r);
  247.       PROPAGATE.left^p(c_plus_p.i);
  248.       IF (DIM1 MOD p) = (DIM1 MOD (2*p)) THEN
  249.         multkomplex(c_plus_p,tmp5,tmp6);
  250.         addkomplex (tmp1,tmp6,c);
  251.         subkomplex (tmp1,tmp6,tmp7);
  252.       END;
  253.       PROPAGATE.right^p(tmp7.r,c_plus_p.r);
  254.       PROPAGATE.right^p(tmp7.i,c_plus_p.i);
  255.       IF (DIM1 MOD p) <> (DIM1 MOD (2*p)) THEN
  256.          c := c_plus_p;
  257.       END;
  258.     ENDPARALLEL;
  259.   END;
  260.   STORE (c,b);
  261.   output (b);
  262. END fastfouriertrans.
  263.  
  264.